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Abstract 

Three-dimensional numerical simulations of decaying turbulence in a magnetized plasma are per- 
formed using a so-called FLR-Landau fluid model which incorporates linear Landau damping and 
finite Larmor radius (FLR) corrections. It is shown that compared to simulations of compressible 
Hall-MHD, linear Landau damping is responsible for significant damping of magnetosonic waves, 
which is consistent with the linear kinetic theory. Compressibility of the fluid and parallel energy 
cascade along the ambient magnetic field are also significantly inhibited when the beta parameter 
is not too small. In contrast with Hall-MHD, the FLR-Landau fluid model can therefore correctly 
describe turbulence in collisionless plasmas such as the solar wind, providing an interpretation for 
its nearly incompressible behavior. 
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I. INTRODUCTION 



Hydrodynamics and Magnetohydrodynamics (MHD) are the central descriptions used to 
study turbulence in the solar wind and in a wide range of natural systems. Specifically for the 
solar wind, MHD description yielded a great success in our understanding of observational 
data (e.g. see reviews by Goldstein et al. [l], Tu and Marsch 2], Bruno and Carbone 



Horbury et al. [4], Marsch 
typical 



Ofman 



Bruno 



. Observational studies show that the solar wind is 
v found to be only weakly compressible (e.g. see Matthaeus et al. |7], Bavassano and 
and the reviews cited above) and usual turbulence models which predict the enerj 



spectra are derived in the framework of an incompressible MHD description (Iroshnikov [9j, 



Kraichnan 



10|, Goldreich and Sridhar (ll|, ll2j, Galtier et al. [13|, |lij, Boldyrev 



dyrev 



19 



Lithwick et al. 17 1, Chandra n [18| , Perez and Bo 
see also reviews by Cho et al. [21], Zhou et al. 22], Galtier 



, Podesta and B 



rattacharjee 



231 ] and Sridhar 



24J). Theoretical 



models which describe the radial evolution of spatially averaged solar wind quantities are also 
usually developed in the framework of incompressible MHD formulated in Elsasser variables 



(e.g. Zhou and Matthaeus 
Matthaeus et al. 



2911 . Smith et al. 



30|, 



25H27|, Marsch and Tu [28|], Zank et al. 
3jJ, Breech et al. {32]). It is however well known that the solar wind 



is not completely incompressible and many phenomena which require compressibility are 
observed in the solar wind, such as the evolution of density fluctuations (e.g. Spangler and 



Armstrong [33J], Armstrong et al. [34j . Coles et al. [35|, Gra 



1 et al. [36j, Woo and Habbal 



40[), magnetic holes, solitons 



Bellamy et al. [38], Wicks et al. [39[, Telloni et al. 
and mirror mode structures (e.g. Winterhalter et al. 411, Franz et al. 421 , Stasiewicz et 

n n 

al. [43j, Stevens and Kasper [44J) or strong temperature anisotropies which trigger, and are 
limited by, micro-instabilities such as mirror and fire-hose (e.g. Gary et al. (45], Kasper 



et al. [46|], Hellinger et al. [47j, Matteini et al. 48], Bale et al. 49]). The importance of 
compressibility was stressed by Carbone et al. J50], who compared the observational data 
with the energy flux scaling laws of Politano and Pouquet 51[, which are exact relations of 
incompressible MHD. Carbone et al. determined that the scaling relations can better fit the 
data if the relations are phenomenologically modified to account for compressibility. The 
incorporation of weakly compressional density fluctuations was partially addressed by so- 
called nearly incompressible models, which expand the compressible equations with respect 
to small sonic Mach number (Matthaeus and Brown 52J, Zank and Matthaeus 53l455j) 
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and which were recently formulated in the presence of a static large-scale inhomogeneous 



background (Hunana and Zank 



561 ] . Hunana et al. |57|, |58|, see also Bhattacharjee et al. 



59|). These models however specifically assume, and cannot explain, why the solar wind is 



only weakly compressible. The theoretical compressible MHD models developed in the wave 
turbulence formalism (e.g. Kuznetsov [60], Chandran [61]) also cannot address this issue. 

Describing the solar wind with fully compressible MHD or Hall-MHD formalisms yields 
several problems. Most importantly, these compressible descriptions introduce sound waves 
and slow magnetosonic waves. As elaborated by Howes [f^J, slow magnetosonic waves are 
strongly damped by Landau damping in the kinetic Maxwell- Vlasov description. The pres- 
ence of fast and slow magnetosonic waves naturally implies higher level of compressibility 
and overestimates the parallel energy transfer. Also, numerical simulations of compressible 
Hall-MHD performed by Servidio et al. jf^ in the context of the magnetopause boundary 
layer showed that compared to the usual turbulence in compressible MHD simulations, which 
consists of Alfven waves, the Hall term is responsible for decoupling of magnetic and velocity 
field fluctuations. In compressible Hall-MHD regime, Servidio et al. observed spontaneous 
generation of magnetosonic waves which transform to a regime of quasi perpendicular "mag- 
netosonic turbulence". This is in contrast with observational studies which typically show 
that turbulence in the solar wind predominantly consists of quasi perpendicular (kinetic) 



Alfven waves (e.g., Bale et al. 64j . Sahraoui et al. 65|, |66|[). Finally, all usual MHD or 
Hall-MHD models cannot address how the energy is actually dissipated at small scales. It 
is well known that solar wind plasma is almost collisionless and therefore the classical von 
Karman picture of energy being dissipated via viscosity is not applicable to the solar wind. 
It is evident that new and more realistic models have to be introduced which can overcome 
some of these drawbacks. 

The most realistic approach is of course the fully kinetic Vlasov-Maxwell description. It 
is however analytically intractable, and even the biggest kinetic simulations cannot resolve 
the large-scale turbulence dynamics. Two leading approaches appear to be promising to 
substitute MHD in describing solar wind turbulence : Gyrokinetics and Landau fluids. 
Gyrokinetics (e.g. Schekochihin et al. Howes et al. {ssj) was originally developed 

for simulations of fusion in tokamaks. It is a kinetic-like description, which averages out 
the gyro-rotation of particles around a mean magnetic field and therefore makes the kinetic 
description more tractable, mostly by eliminating fast time scales. Derived directly from the 
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kinetic theory, gyrokinetics has a crucial advantage of being asymptotically correct. Landau 
fluid description, on the other hand, is a fluid-like extension of compressible Hall-MHD, in 
which wave dissipation is incorporated kinetically by the modeling of linear Landau damping, 
thus retaining a realistic sink of energy. Other linear kinetic effects such as finite Larmor 
radius corrections are also incorporated. 

The simplest Landau fluid closure was considered by Hammett and Perkins jf^]] and the 
associated dispersion relations were numerically explored by Jayanti, Goldstein and Vinas 



TCH ] . The Landau fluid model was further advanced by Snyder, Hammett and Dorland 71] . 



who considered the largest MHD scales, starting from the guiding center kinetic equation. 
The approach was reconsidered and refined to its present form with incorporation of Hall 
;erm and finite Larmor radius (FLR) corrections by Passot, Sulem, Goswami and Bugnon 



72 



475|. This Landau fluid approach starts with the Vlasov-Maxwell equations and derives 
nonlinear evolution equations for density, velocity and gyrotropic pressures. In the simplest 
formulation the model is closed at the level of heat fluxes by matching with the linear 
kinetic theory in the low frequency limit. Kinetic expressions usually contain the plasma 
dispersion function which is not suitable for fluid-like simulations. Landau fluid closure is 
therefore performed in a way as to minimize occurrences of this plasma dispersion function 
and, where not possible, this function is replaced by a Pade approximant. This eliminates 
the time non-locality and also results in the presence of a Hilbert transform with respect 
to the longitudinal coordinate (in the direction of the ambient magnetic field) which in 
the fluid formalism is associated with linear Landau damping. Further details about the 
development of Landau fluid models are thoroughly discussed in the papers cited above. 
The Landau fluid approach should however be contrasted with the more classical gyrofluid 
models (e.g. Dorland and Hammett }7^ |. Brizard {7?]], Scott which are derived by 

taking fluid moments of the gyrokinetic equation and for which a similar closure scheme is 
applied afterwards. 

The Landau fluid approach has the following advantages. In contrast with Hall-MHD, 
it contains separate equations for parallel and perpendicular pressures and heat fluxes. It 
therefore allows for the development of temperature anisotropy, which is observed in the 
solar wind. Noticeably, in contrast with gyrokinetics or gyrofluid models, the Landau fluid 
model does not average out the fast waves. Compared with these approaches, it also has 
an advantage in that the final equations including the FLR corrections are written for the 
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usual quantities measured in the laboratory frame. Existing spectral MHD and Hall-MHD 
codes can therefore be modified to Landau fluid description relatively easily. Also impor- 
tantly, even though gyrokinetics is a reduced kinetic description, it is still 5-dimensional and 
therefore naturally quite difficult to compute. While current largest numerical simulations 
of gyrokinetics (Howes et al. 79|) require thousands CPU-cores for a fluid-like 128 x 64 2 res- 
olution, the FLR-Landau fluid model requires computational power only slightly larger than 
the usual Hall-MHD simulations. The results presented here, which employ a resolution of 
N = 128 grid points in all three directions, were calculated using 32 CPU-cores. 

Landau fluid models can be developed with several levels of complexity. For these first 
Landau fluid simulations of three dimensional turbulence, we use a simplified version of the 



most general Landau fluid model 



72], where we constrain ourselves to isothermal electrons 



and leading order corrections in terms of the ratio of the ion Larmor radius to the considered 
scales. A similar model was used by Borgogno et al. [8fJ, who studied the dynamics 
of parallel propagating Alfven waves in a medium with an inhomogeneous density profile. 
They numerically showed that the observed Alfven wave filamentation and later transition 
to the regime of dispersive phase mixing is consistent with particle-in-cell simulations. In 
this paper we concentrate on freely decaying turbulence and compare these to simulations 



of compressible Hall-MHD. Numerical integration o: 



the full Landau fluid model in one 



8jJ , who investigated the dynamics very 



space dimension was presented by Borgogno et al. 
close to the mirror instability threshold and showed the presence of magnetic holes. Results 
considering quasi-transverse one-dimensional propagation in the full Landau fluid model are 



presented in [82 



II. THE MODEL AND ITS NUMERICAL IMPLEMENTATION 

Considering a neutral bi-fluid consisting of protons (ions) and isothermal electrons, the 
Landau fluid model consists of evolution equations for proton density p = m p n (where 
m p is the proton mass and n the number density), proton velocity u p , proton parallel 
and perpendicular pressures p\\ p , p± p and heat fluxes g^, q± p , together with the induction 
equation for magnetic field b. The equations are normalized and density, magnetic field 
and proton velocity are measured in units of equilibrium density po, ambient magnetic field 
B , and Alfven speed Va = B / y/47ip , respectively. Pressures are measured in units of 
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initial proton parallel pressure p^ and heat fluxes in units of p\\j}Va- The total pressure 
in the momentum equation has the form of a tensor. Defining b = b/\b\ as a unit vector 
in the direction of local magnetic field, the proton pressure can be cast in the form p p = 
pj_ p n + p\\ p T + II, where r = b <g> b, and n = I — b <g> b, with i" being the unit tensor. Finite 
Larmor radius corrections to the gyrotropic pressures are represented by II. Operator ® 
represents the usual tensor product and in the index notation, for example, Tjj = bibj. 
Electrons are assumed to be isothermal with the scalar pressure p e = nT e (0) , where T e (0) is 
the electron temperature. Parameter R iy whose inverse multiplies the Hall-term and also 
the FLR corrections, is defined as Ri = L/di, where di is the ion inertial length and L is 
the unit length. The proton plasma beta is defined with respect to parallel pressure and 
= 8ttp^/Bq. The density, momentum and induction equations of the FLR-Landau fluid 
model can then be expressed as 

^ + V • (pu p ) = 0, (1) 
+ u p ■ Vu p + ^ V • (p ±p n + p\\ p r + n + p e I) 

--(V x b) x b = 0, (2) 

P 



- = VxKxb)--Vx 



-(V x b) x b 

P 



(3) 



Dropping, for simplicity, indices p for proton velocity u p and proton pressures p± p ,p\\ p , the 
evolution equations for perpendicular and parallel pressures reads (neglecting the work done 
by the FLR stress forces) 

^ + V • (p ± u) + P± V u-p ± b- (Vu) • b 

+V • (q ± b) + q ± V-b = 0, (4) 
^ + V • + 2 P ||b • (Vu) ■ b + V • (q\\b) - 2q ± V -6 = 0. (5) 

Assuming an ambient magnetic field of amplitude B in the positive z-direction, a semi- 
linear description of the finite Larmor radius corrections in the pressure tensor neglecting 
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heat flux contributions can be expressed as 

<P±> 



2R { 

(P±> 



2-R,; 



U yz = U zy = — [2{p\\)d z u x + (p±)(d x u z - d z 

its. 



u, 



n xz = U zx = -— [2(p\\)d z u y + (p±){d y u z - d. 



n Z2 = o, 



(6) 

(7) 

(8) 

(9) 
(10) 



where (p±) and (p\\) represents the instantaneous averaged ion pressures over the entire 
domain, whose time variation is aimed to take into account the evolution of the global prop- 
erties of the plasma. Finally, parallel and perpendicular heat fluxes q\\, q± evolve according 
to 




1 g| 
!-¥ 2 



dz(P\\~ P); 



(11) 



h(0) 

lb 



(0) 



(12) 



where rlpjTiip- 1 are the initial perpendicular and parallel proton temperatures and d/dt is 
the convective derivative. The operator %, which is defined as 

30 f(z') 



Hf{z) 



1 



VP 



-dz'. 



z — z 



(13) 



reduces in the Fourier space to a simple multiplication by ik z /\k z \ and, is the signature of 
the linear Landau damping. 

To gain physical insight into a quite complicated model ([!])- (TP2"]) . it is useful to momentar- 
ily consider just the largest scales by putting 1/ Ri — >■ 0, which eliminates the nongyrotropic 
contributions to the pressure tensor and which also eliminates the Hall term. The result- 
ing set of equations still contains the linear Landau damping and with the exception of 
isothermal electrons, it is analogous to the model of Snyder, Hammett and Dorland 7l[. If 
this model is further simplified by elimination of eq. ( TTT]) . (fT2j) and by instead prescribing 
1\\ — Q± — 0, the resulting model collapses to the double adiabatic model (Chew et al. [84j , 



see also Kulsrud 



851 ] ) and the Landau damping disappears. The presence of ion Landau 



damping in the system (fT|)- (IT2j) is therefore a result of closure equations (ITT)) . (jT2j) for the 
heat fluxes, which contain the operator TL. At least in the static limit, this closure can be 
viewed as a modified Fick's law where the gradient operator that usually relate the heat flux 
to the temperature fluctuations is here replaced by a Hilbert transform that is a reminis- 
cence of the plasma dispersion function arising in the linear kinetic theory, and is a signature 
of Landau (zero-frequency) wave-particle resonance. The effect of Landau damping in the 
system (TOl- fTS!) might be better understood by solving the linearized set of equations. In 
the Appendix we consider waves which propagate parallel to the ambient magnetic field. 
It is shown that except the Alfven waves, linear waves have a frequency with a negative 
imaginary part, and are therefore damped. This corresponds to linear Landau damping. 

The model used for the simulations presented here has several limitations. First of all, it 
contains electrons which are assumed isothermal, a regime in fact often assumed in hybrid 
simulations which provide a kinetic description of the ions and a fluid description of the 
electrons. In the future, more realistic simulations will be performed with inclusion of 
independent evolution equations for parallel and perpendicular electron pressures and heat 
fluxes. This will also result in the presence of electron Landau damping, which is absent in 
the model presented here and which seems to play an important role in solar wind turbulence. 
Another main limitation appears to be the form of finite Larmor radius corrections, which 
are derived as a large-scale limit of FLR corrections of the full Landau fluid model and 
which are therefore significantly simplified. The FLR corrections (I6l)-(fl0l) are sufficient for 
the simulations of freely-decaying turbulence, which do not lead to significant temperature 
anisotropies. However, our preliminary simulations which employ forcing and lead to strong 
temperature anisotropies show that the FLR corrections (l6l)- (IT0l) are overly simplified and 
might lead to artificial numerical instabilities. For simulations with strong temperature 
anisotropies, a more refined description of FLR corrections is required (see Passot and Sulem 
721 ] . Borgogno et al. [81]). 

To explore the behavior of the FLR-Landau fluid model, we performed simulations of 
freely decaying turbulence. The code we used is based on a pseudo-spectral discretization 
method, where spatial derivatives are evaluated in Fourier space. The time stepping is 
performed in real space with a 3rd order Runge-Kutta scheme. Spatial resolution is N 3 = 
128 3 and the size of the simulation domain is L = 16 x (2tt) in each direction. The Hall 



S 



parameter is Ri = 1, implying that the lengths are measured in the units of ion inertial 
length di. Velocity and magnetic field fluctuations of mean square root amplitude (u 2 ) 1 / 2 = 
(6 2 ) 1 / 2 = 1/8 are initialized in Fourier space on the first four modes kdi = m/16, where 
m G [1,4], with flat spectra and random phases and are constructed to be divergence free. 
Constant pressures p± = p\\ = 1, density p = 1 and heat fluxes q± = q\\ = are initialized 
in the entire domain. The temperature of the electrons is = 1 and is therefore equal 
to proton temperatures, which can be defined as T± = (p±/p) and Tji = (p\\/p)- In the 
simulations presented here, both T± and Tji stay rather close to their initial value. The 
compressible Hall-MHD model consists of eq. ([I])-®, where the divergence of the pressure 
tensor in eq. (J2]) is substituted with the usual gradient of scalar pressure and, assuming 
the adiabatic law p = p 7 , in the normalized units it is equal to /3 /7Vp 7 , where (3 is the 
usual plasma beta defined as (3 = c 2 /Vj and the sound speed c 2 = jp/p. Adiabatic index 
7 = 1.66 is used. 



As shown for example by Hirose et al. 86[ and Howes 62|, it is not straightforward to 



compare Hall-MHD and Vlasov-Maxwell kinetic theory because the associated dispersion 
relations are quite different if in the kinetic description the proton (ion) temperature T p 
is not negligible with respect to the electron temperature T e . We choose to follow Howes 



621 ] . who, in order to compare Hall-MHD with the kinetic theory (which is represented 



here by the FLR-Landau fluid model), defined the necessary relation between /3 and /3n 
as /3q = /3||(1 + Tp/T e )/(2T P /T e ). For equal proton and electron temperatures T p = T e this 
relation yields (3o = (3\\. Landau damping alone is not sufficient to run the code and some 
kind of artificial dissipation is needed to terminate the cascade. We here resorted to use a 
filtering in the form of |{1 — tanh[(m — 0.8iV/4)/3]}, where m is the mode index, applied 
each time step on all fields for both Hall-MHD and Landau fluid regimes. Because of the 
filtering, it is crucial to have identical time steps dt = 0.128 in both models. In most of 
the simulations we used {3o = (3\\ = 0.8. In normalized units, the Alfven speed is equal to 
unity, the usual MHD sound speed c s = a/5o = 0.894 and the turbulent sonic Mach number 
M s = (it 2 ) 1 / 2 jc s = 0.14. In section IV., we also consider simulations with (3 = j3\\ = 0.25 
(M B = 0.25) and (3 = (3» = 0.1 (M s = 0.4). 
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frequency (co) frequency (ra) 



FIG. 1: Left and right polarized Alfven waves for the propagation angle of 0° for Hall-MHD 
(left) and Landau fluid (right), as revealed by frequency analysis of b x modes with wavenumbers 
k x = 0, k y = 0,k z di = m/16, where m = 1 (red), m = 2 (green), m = 4 (blue), m = 8 (black). 
Theoretical predictions for peaks calculated from dispersion relations for given k are shown on the 
top axis. 

III. IDENTIFICATION OF THE MHD MODES 

A wave analysis procedure was implemented in the code, which consists in choosing few 
modes in spatial Fourier space for each field and recording their value every 20 time steps. 
After the run, the time Fourier transform is performed and frequency-power spectra are 
obtained for each mode. This procedure makes it possible to identify at which frequency 
there is maximum power for each mode, and by comparing these with frequencies obtained 
from theoretical dispersion relations it allows to uniquely identify which waves are present 
in the system. This method was previously used for incompressible MHD by Dmitruk and 



Matthaeus 87], therefore detecting Alfven waves. Considering the dynamics associated 
with waves propagating in the direction parallel to the ambient magnetic field (propagation 
angle 6 = 0°), Fig. 1 shows frequency-power spectra of four modes with wavenumbers 
k x = 0, k y = 0, k z d{ = m/16, with m = 1,2, 4, 8, recorded from the component b x for Hall- 
MHD (left) and for FLR-Landau fluid (right). For Hall-MHD, these correspond to polarized 
Alfven waves which obey the dispersion relation u = ±k 2 /(2Ri) + ky/l + (k/2Ri) 2 . The 
theoretical frequency values which are expected from this dispersion relation for given k are 
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plotted on the top axis of the figure as black vertical lines. Figure 1 shows that the resolution 
of 128 3 is sufficient to clearly distinguish between left and right polarized Alfven waves. In 
contrast with Hall-MHD, for which it is possible to write the general dispersion relation 
for frequency u as a relatively simple polynomial of 6th order, for FLR-Landau fluids the 
general dispersion relation would be uneconomically large to write down. In general, it is 
necessary to numerically solve the determinant obtained from linearized equations (IXT)- fTT2|) 
for a given wavenumber k after assuming linear waves. FLR-Landau fluid ffTj)- ffT2]) consists 
of 11 evolution equations in 11 variables and, together with the divergence free constraint 
for the magnetic field, therefore yields general dispersion relation for frequency u in the 
form of a complex polynomial of 10th order. This represents 5 forward and 5 backward 
propagating waves, with some solutions having highly negative imaginary part and which 
are therefore strongly damped. A similar situation is encountered in the Vlasov-Maxwell 
kinetic theory which essentially yields an infinite number of strongly damped solutions. For 
the propagation angle 8 = 0° and additional constraint T± = Tji = 1, it is however possible 
to obtain an analytic solution for the circularly polarized Alfven waves as 

with two other solutions obtained by substituting oj with —oj. The dispersion relation is 
quite similar to that of Hall-MHD with additional terms proportional to j3\\ and resulting 
from the finite Larmor radius corrections. Expected theoretical frequencies obtained from 
this analytic solution are plotted on the top axis of Fig. 1 for the Landau fluid regime 
(right). They again match quite precisely. Note also the moderately strong damping of 
parallel Alfven waves in Landau fluid regime, which can be seen for the last mode m — 8. 
Landau damping does not act directly on linear constant amplitude Alfven waves obeying 
relation (Tl4"|) . which are exact solutions of linearized FLR-Landau fluid model. However, 
nonlinear parallel Alfven waves, which are of course present in the full model, cause produc- 
tion of density (sound) fluctuations. Sound waves in the FLR-Landau fluid model, as well 
as in the kinetic theory, are heavily damped by Landau damping as shown below and this 



process therefore also results in damping of Alfven waves. Mj0lhus and Wyller [88] stud- 
ied the kinetic derivative nonlinear Schrodinger equation (KDLNS) for parallel propagating 
long-wavelength Alfven waves where they refer to this effect as nonlinear Landau damping, 
because it is acting on Alfven waves in a nonlinear way. 
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FIG. 2: Sound waves for the propagation angle of 0° for Hall-MHD (left) and Landau fluid (right), 
from frequency analysis of u z modes with the same wavenumbers as in Fig. 1. Sound waves (S) 
are heavily damped for the Landau fluid, which is consistent with the kinetic theory. The spectra 
also show weak presence of Alfven waves (A), which are visible for the first two modes. 

Further exploring propagation angle of 9 = 0°, Fig. 2 shows frequency power spectra 
obtained from component u z and which should therefore predominantly display sound waves 
(almost identical spectra can be obtained from component p). The same modes with k z di = 
m/16, where m = 1,2,4,8 are shown as in Fig. 1. For Hall-MHD, parallel sound waves 
obey the dispersion relation oj = kc s , where the sound speed c s = vT^o- Sound waves (S) 
are clearly presented in Fig. 2 for Hall-MHD (left) with quite sharp peaks which match 
the theoretical dispersion values shown on the top axis. Weak presence of polarized Alfven 
waves (A) for modes m — 1,2 is also visible in this component generated by nonlinear 
coupling. For the Landau fluid regime (Fig. 2 right), it is not possible to obtain simple 
analytic dispersion relation which corresponds to sound waves and correct values must be 
obtained numerically as explained above (see also the Appendix). This yields 5 frequencies 
for forward propagating waves, with 2 solutions corresponding to the polarized Alfven waves 
(I14p and 3 solutions which are highly damped. The sound wave was chosen as the least 
damped solution. Dispersion relation uj = ky~(3~+ )/3\\/2 obtained from the double 
adiabatic model, to which Landau fluid description collapses after assumption of zero heat 
fluxes and large scales, can also be used as a heuristic guide to determine which out of 
the 3 damped solutions represents the sound wave frequency. Obtained frequency values 
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are again plotted on the top axis. Figure 2 shows that the sound waves (S) are heavily 
damped for the Landau fluid regime (Fig. 2 right) with the nonlinearly generated Alfven 
waves (A) completely overpowering the spectra. Inhibition of sound waves is consistent with 



the kinetic theory. As elaborated by Howes 



62], sound waves are overestimated by MHD 



and Hall-MHD descriptions as they represent solutions which are strongly damped by the 
Landau resonance. Note that the calculated sound wave frequencies in the FLR-Landau 
fluid regime are actually higher than the corresponding Alfven wave frequencies. A similar 
result was obtained by Howes 62| who numerically compared dispersion relations for Hall- 
MHD and kinetic theory and who, for the nearly parallel propagation with Tj = T e and 
0o — P\\ = 1, noted that the kinetic solution corresponding to the slow wave has a higher 
phase speed than the kinetic solution corresponding to the Alfven wave. 

Considering the propagation angle of 9 = 45°, Fig. 3 shows frequency-power spectra in 
component b x for wavenumbers k y = 0, k x di = k z di = m/16, where m = 1, 2, 4, 6 and which 
therefore predominantly shows slow (S) and fast (F) magnetosonic waves. Spectra also show 
weaker presence of Alfven waves (A). Again for this angle of propagation, the slow waves 
are strongly damped. In both Hall-MHD and FLR-Landau fluid simulations, the associated 
dispersion relations had to be solved numerically and the predicted frequencies for spectral 
peaks are shown on the top axis. Slow waves (S) in the FLR-Landau fluid regime were again 
identified as the least damped frequency out of the 3 highly damped solutions. For identical 
wavenumbers, the presence of Alfven waves can be better explored in the component b y and 
corresponding spectra are shown in Fig. 4. 

Finally, considering purely perpendicular propagation with 9 = 90°, Fig. 5 shows 
frequency-power spectra obtained from the density field p and wavenumbers k y = 0, k z = 
0,k x di = m/16, where m = 1,2,4,8. For Hall-MHD regime (Fig. 5 left) the spectral peaks 



correspond to magnetosonic waves with the usual dispersion relation u = kyl + /3q. In 
the Landau fluid regime, linearized set of equations (fTT)-(IT2l) with the additional constraint 
T± — Tji = 1 can be shown to yield the dispersion relation for perpendicular magnetosonic 
waves in the form 



bj = k 



\ 



The dispersion relation clearly shows the effect of inclusion of isothermal electrons (for 
simulations presented here Tj ^ = 1) and also the effect of finite Larmor radius corrections 
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FIG. 3: Slow (S) and fast (F) magnetosonic waves for the propagation angle of 45° for Hall- 
MHD (left) and Landau fluid (right), from frequency analysis of b x modes with wavenumbers 
k y = 0, k x di = k z di = m/16, where m = 1 (red), m = 2 (green), m = 4 (blue), m = 6 (black). The 
presence of Alfven waves (A) is also visible. Theoretical predictions for slow and fast waves are 
shown on the top axis. 

which are represented by the last quadratic term. Theoretical predictions from Hall-MHD 
and FLR-Landau fluid dispersion relations are again shown on the top axis. To clearly show 
the shift of the peaks between the two regimes, we also added the theoretical Hall-MHD 
frequencies to the top axis of FLR-Landau fluid regime (Fig. 5 right) and represent them 
with the small magenta lines. 

IV. FLOW COMPRESSIBILITY 

Compressibility of the flow can be evaluated by decomposing the velocity field into its 
solenoidal and non-solenoidal components and by calculating the associated energies accord- 
ing to 

\ Uk \ =2^ | fc i2 +!L puis ' ( 16 ) 

k k k 

where the left-hand side corresponds to the total energy E u in velocity field, the first term in 
the right-hand side corresponds to the energy Ei n in the solenoidal component and the second 
term E c originates from the compressible one. Relation (fl~6j) can be therefore expressed as 
E u = E in + E c and the compressibility of the flow can be evaluated as a ratio of compressible 
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FIG. 4: Alfven waves for the propagation angle of 45° for Hall-MHD (left) and Landau fluid 
(right), from frequency analysis of b y modes with the same wavenumbers as in Fig. 3. Theoretical 
predictions for the Alfven waves frequencies are shown on the top axis. 

and total energy E c /E u . Time evolution of EJE U for Hall-MHD and FLR-Landau fluid 
regime with /3q = (3\\ = 0.8 is presented in Fig. 6. The figure shows that the ratio E c /E u 
which represents compressibility is significantly lower in Landau fluid regime and is therefore 
a result of presence of Landau damping. The question then arises of the influence of the 
sonic Mach number in the compressibility evolution. The time evolution of E c /E u for 
A) = 0\\ — 0.25 (which corresponds to M s = 0.25) is displayed in Fig. 7 (left), whereas 
the time evolution of E c /E u for (3 = /3n = 0.1 (which corresponds to M s = 0.40) is shown 
in Fig. 7 (right). The simulations were performed with the same time step dt = 0.128 
and filtering, for approximately half the total integration time of the previous regime with 
0o — P\\ — 0.8 (M s = 0.14). The figure shows that while the compressibility is still clearly 
reduced in Landau fluid simulation with /3 = 0\\ = 0.25, this reduction is almost insignificant 
for simulations with fto = (3u = 0.1. This is an expected effect, as the strength of the Landau 
damping is proportional to /V We note that the turbulent sonic Mach number in the 
solar wind is typically small and, for example, analysis of observational data performed by 
Bavassano and Bruno js] (from 0.3-1.0 AU) showed that the most probable value is between 
M s = 0.1 — 0.2 with the distribution having an extended tail to lower values. 

The question also arises of the compared evolution of total energy for Hall-MHD and 
FLR-Landau fluid simulations. The total energy E tot can be evaluated as the sum of the 



15 




10" 

frequency (a>) 




10" 

frequency (ra) 



FIG. 5: Magnetosonic waves for the propagation angle of 90° for Hall-MHD (left) and Landau fluid 
(right), from frequency analysis of density modes with wavenumbers k y = 0, k z = 0, k x di = m/16, 
where m = 1 (red), m = 2 (green), m = 4 (blue), m = 8 (black). Theoretical predictions from 
dispersion relations are shown on the top axis (long black lines). For comparison, on the right 
panel we also included the frequency predictions from Hall-MHD (small magenta lines). 




FIG. 6: Compressibility for Hall-MHD (red line) and FLR-Landau fluid (blue line) evaluated as 
(Sfc \k ' u k\ 2 /\k\ 2 )/ J2k \ u k\ 2 for /3rj = /3|| = 0.8. Both regimes start with the identical initial 
condition where the velocity field is divergence free. The figure shows that the compressibility is 
clearly inhibited in the Landau fluid simulation. 
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FIG. 7: Compressibility for Hall-MHD (red line) and FLR-Landau fluid (blue line) evaluated as 
(£ fc |fc • Wfc| 2 /I fc | 2 )/Efc \ u k\ 2 with /3 = /3|| = 0.25 (left) and with /3 = P\\ = 0.1 (right). The 
figure shows that the compressibility is clearly inhibited in Landau fluid for simulations with 
A) = P\\ = 0.25, whereas for simulations with (3q = /3| =0.1 the level of compressibility is almost 
identical. 

kinetic energy Eki n , the magnetic energy E mag and the internal energy Ei nt . In Hall-MHD 
and FLR-Landau fluid model, the definitions of kinetic and magnetic energy are identical 
and equal to 

E k i n = ^ J p\u\ 2 dx 3 , E mag = ~J \b\ 2 dx 3 . (17) 

However, the definition of internal energy is naturally different in each model. In the Hall- 
MHD model, the internal energy is defined as 

HMHD: E mt = - j P Ux\ (18) 

whereas in the FLR-Landau fluid model with isothermal electrons, the internal energy is 
given by 

Landau: E mt = 4 J (p ± + ^ + r^plnp) dx\ (19) 

It is emphasized that because of the filtering, the total energy E tot is not exactly preserved 
in the Hall-MHD and Landau fluid simulations. During the simulations with (3$ = (3\\ = 0.8, 
the total energy in Hall-MHD decreased by approximately 1.2%, whereas in FLR-Landau 
fluid the decrease was approximately 0.5%. Filtering indeed dissipates kinetic and magnetic 
energies without turning them into heat. Therefore, in contrast with Landau fluid models, 
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FIG. 8: Normalized mechanical fluctuations energy E^ n + E mag (left) and normalized internal 
energy Ei n t (right) for Hall-MHD (red line) and FLR-Landau fluid (blue line), in the case /3n = 
/9ii = 0.8. Landau damping transfers energy from E^ + E mag into Ei nt . 

Hall-MHD simulations do not have heating at all and the internal energy could increase 
only through development of density fluctuations. Time evolution of the sum of kinetic and 
magnetic energies (normalized to their initial values) is displayed in Fig. 8 (left), where the 
energy contained in the ambient magnetic field was subtracted. This figure shows that this 
sum decays faster for FLR-Landau fluid simulation. Time evolution of internal energy Ei nt 
is shown in Fig. 8 (right), where energies were again normalized to their initial value. The 
initial jump observed in both simulations reflects a rapid adjustment from initial conditions 
that are not close to an equilibrium state. Later on, the internal energy of Hall-MHD 
simulation decreases, whereas the internal energy of FLR-Landau fluid simulation more or 
less smoothly increases until around time t = 1500. During this time (which corresponds 
to over 10 4 time steps) the Landau damping acts strongly and, by mainly damping slow 
waves, converts the mechanical energy E^ in + E mag into the internal energy E int , which 
represents heating of the plasma. The question also arises what fraction of mechanical energy 
is dissipated directly by Landau damping and what fraction is dissipated by the filtering 
process. Unfortunately, we are unaware of any technique how to address this question. 

We note that for simulations of freely decaying turbulence the heating is quite weak, im- 
plying that driving the system is necessary to produce significant temperature anisotropies. 
However, the absence of forcing, which yields only a small amount of heating, makes it easier 
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to precisely identify various waves in the system as was presented in Sec. III. 

V. ANISOTROPY OF THE ENERGY TRANSFER 

The presence of Landau damping can also be seen in the usual wavenumber velocity and 
magnetic field spectra. Considering first simulations with (3 = (3\\ = 0.8, Fig. 9 shows the 
velocity spectra with respect to perpendicular and parallel wavenumbers k±, k\\ which are 
defined as E u = J E u (k±)dk^ = J E u (k\\)dk\\. With respect to k±, the spectra for Hall- 
MHD and FLR-Landau fluid are almost identical (Fig. 9 left) whereas with respect to kn 
(Fig. 9 right), the spectra of Landau fluid are much steeper. Landau damping therefore 
significantly inhibits the parallel transfer. Even though low resolution does not allow to 
precisely identify the slopes of the spectra, three straight lines were added to figures and 
correspond to power law solutions k s , where s = —3/2,-5/3 and —7/3. For E u (k^_), the 
closest spectral index value appears to be —5/3, the spectral range being however quite 
limited. The same conclusion with the inhibition of parallel transfer is also obtained for 
the magnetic field spectra, which are almost identical to the velocity spectra and are shown 
in Fig. 10. In contrast, for simulations with /3q = P\\ =0.1 the parallel spectrum E(k\\) 
displays a similar behavior for Hall-MHD and Landau fluids. The velocity spectra are shown 
in Fig. 11 and the magnetic field spectra are shown in Fig. 12. This is consistent with results 
presented in the previous section where it was shown that the Landau damping is responsible 
for significant reduction of compressibility for simulations with (3q = (3\\ = 0.8, whereas for 
simulations with (3 = (3\\ = 0.1, the Landau damping was much weaker and the reduction 
of compressibility almost negligible. 

It is useful to compare our compressible simulations to incompressible simulations. Nat- 
urally, our compressible Hall-MHD code cannot be run in an incompressible regime, for 
which the turbulent sonic Mach number M s — > and the sound speed c s — > oo. Neverthe- 
less, incompressible MHD simulations of decaying turbulence were performed, for example, 
by Bigot et al. js^]. These simulations showed that the combined velocity and magnetic 
field spectra (they used Elsasser variable z + ) are much steeper with respect to k\\ than with 
respect to k±, if the ambient magnetic field is sufficiently strong. Our compressible sim- 
ulations presented here have initially (b 2 ) 1 ^ 2 /B = 1/8 and can therefore be compared to 
incompressible MHD simulations of Bigot et al. [89| (their Fig. 17 c,d), where this ratio 
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FIG. 9: Velocity spectra for Hall-MHD (red) and Landau fluid (blue) with respect to perpendicular 
wavenumber E u (k±) (left) and with respect to parallel wavenumber E u (k^) (right), for (3q = = 
0.8. Spectra were taken at time t = 5248. Straight lines correspond to /c -3 / 2 , A; -5 / 3 and A; -7 / 3 . The 
figure shows that spectra with respect to fen are much steeper in Landau fluid simulation, which is 
a result of Landau damping. 



is 1/5 and 1/15, respectively. Considering compressible Hall-MHD, our simulations showed 
that spectra with respect to fcii (e.g. Fig. 9 right, red line) are steeper than spectra with 
respect to k± (Fig. 9 left, red line). However, these parallel spectra are nowhere near as 
steep as the parallel spectra of Bigot et al. [89[. Interestingly, their parallel spectra more 
resemble our parallel spectra for Landau fluid model (Fig. 9 right, blue line), where the Lan- 
dau damping was strong = 0.8). These results show that magnetosonic waves (which 
are not present in an incompressible MHD description) play an important role in regulating 
the parallel energy cascade. Steeper parallel spectra in Landau fluid model are therefore a 
consequence of damping of slow magnetosonic waves by Landau damping. 



VI. CONCLUSION 



We have presented the first three-dimensional fluid simulations of decaying turbulence 
in a collisionless plasma in conditions close to the solar wind. For this purpose, we used 
the FLR-Landau fluid model that extends compressible Hall-MHD by incorporating low- 
frequency kinetic effects such as Landau damping and finite Larmor radius corrections. It 
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FIG. 10: Magnetic field spectra for Hall-MHD (red) and Landau fluid (blue), for E b (k ± ) (left) and 
E h {k\\) (right) when /3 = /3y = 0.8. Spect ra were taken at time t = 5248. 




FIG. 11: Velocity spectra for Hall-MHD (red) and Landau fluid (blue) with respect to perpendicular 
wavenumber E u (k±) (left) and with respect to parallel wavenumber E u (kn) (right) for /3o = /3ii = 
0.1. Spectra were taken at time t = 5248. Straight lines correspond to k~ 3 ^ 2 ,k~ 5 ^ 3 and k~ 7 / 3 . 

was shown that in spite of the turbulent regime, it is possible to precisely identify linear waves 
present in the system. Comparisons between compressible Hall-MHD and FLR-Landau 
fluid model showed that when beta is not too small, linear Landau damping yields strong 
damping of slow magnetosonic waves in Landau fluid simulations. These waves are indeed 
damped in kinetic theory described by the Vlasov-Maxwell equations but not in compressible 
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FIG. 12: Magnetic field spectra for Hall-MHD (red) and Landau fluid (blue), for E b (k ± ) (left) and 
E b (k\\) (right) when (3 = /3\\ = 0.1. Spectra were taken at time t = 5248. Spectra for Hall-MHD 
and Landau fluid model are again almost identical. 

MHD and Hall-MHD descriptions, which overestimate compressibility and parallel transfer 
in modeling weakly collisional plasmas. The FLR-Landau fluid model can therefore be useful 
for simulating the solar wind, which is typically found to be only weakly compressible. 
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APPENDIX 

To clearly understand how the Landau damping acts in the present system (1)-(12), it is 
useful to solve dispersion relations for linear waves propagating in parallel direction to the 
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ambient magnetic field. A detailed analysis of linear waves for various propagation angles 
was elaborated by Passot and Sulem 



proton temperature anisotropy as T^/T^ = a p , and the normalized electron temperature 



90( . To simplify the analytic expressions, we define the 

(0) 

\p 

as /Ty = t. It can be shown that for parallel propagation angle, the Landau fluid 
model contains four dispersive Alfven waves. Two waves obey the dispersion relation which 
can be expressed as 



k 

LO = 



2Ri 



1 + +^/l + ?K-l)+(^J 1 — ( 1 — ^ ) , (20) 




with another two solutions obtained by substituting u with —u. Obviously, these Alfven 
waves are independent of the electron temperature r, which is a consequence of electrons 
being modeled as isothermal. For a more general Landau fluid model which contains evo- 
lution equation for electron pressures and heat fluxes, the electron temperature r enters 
the dispersion relation for Alfven waves. The solutions (120]) can become imaginary, if the 
expression under the square root becomes negative. At large scales (when 1/Ri — > 0) the 
condition 1 + fi\\{a p — l)/2 < represents the well known criterion for fire-hose instability 
(see, for example, Ferriere and Andre 911]). The Hall term and FLR corrections modify the 
instability criterion. For isotropic temperatures (a p = 1), the solution ( 120]) naturally col- 
lapses to the solution (THl) . The four Alfven waves (1201) can be eliminated from the general 
dispersion relation and this yields a complex polynomial of 6th order in frequency u. Solu- 
tions of this polynomial represent 3 forward and 3 backward propagating waves which have 
a negative imaginary part and are therefore damped. Importantly, it is possible to eliminate 
the dependence on (3\\ and wavenumber k and, after applying a substitution Q = uj/(ky/]3\\), 
the polynomial of 6th order can be simplified to 

0> + n y?- 4 + 3 ; + ii .-^4 + 3 : )+24 + 8r + ^ r - .(9 + 3r)/4 



-16 + 67T -16 + 6tt v -8 + 3vr 

,(15 + 5r)/4 - 2 (5 + 3r)/2 _ = 

-8 + 3tt v -8 + 3tt -8 + 3tt v ; 

Because this polynomial in Q does not depend on j3\\ or k, the substitution implies that 
all 6 waves are linear with k and \J~fi\\- The polynomial ( )2T|) has to be solved numeri- 
cally for a given value of r. The simulations presented here use r = 1 and numerically 
solving polynomial (EQ yields Vt = ±1.48106 - z0.36117, ft = ±0.65467 - i0.88285 and 
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n = ±0.55104 — iO. 44311. These six waves therefore satisfy 

u = fcy^j"(±1.48106-i0.36117), (22) 
u = £^(±0.65467 - i0.88285), (23) 
u = A;^(±0.55104-i0.44311). (24) 

The least damped solution ( 122]) represents the sound wave. The solutions of Landau fluid 
model (1)-(12) for parallel propagation angle are therefore 4 Alfven waves (1201) . 2 sound 
waves ( 122]) and 4 waves ( 123]) . (1241) . which are highly damped. These 4 waves (123]) . (124]) 
do not have an analogy in Hall-MHD description and correspond to solutions of kinetic 
Maxwell- Vlasov description, which contains an infinite number of highly damped solutions. 
Interestingly, the last solution (124]) is not dependent on the value of r and it can be expressed 



analytically as Q = ±\/8 — 7r/4 — 2a/7t/4. After eliminating these waves from eq. (121"]) . the 
polynomial which contains the sound waves (122]) and solutions (123}) is now of 4th order in Q 
and expressed as 

n . + < _V!_ n ,_19.-16 + (-8 + a,)r tf _ ; v^(3 + r) „ + gjl + r) 

-8 + 3tt 2 -8 + 3tt -8 + 3tt -8 + 3tt v ; 

It is of course possible to use Ferrari-Cardano's relations to solve this polynomial analytically, 
the final result is however too complicated and it is still more convenient to solve (125]) 
numerically for a given value of r. 

If this wave analysis is repeated with the model (1)-(10) with heat flux equations q\\ = 0, 
q± = 0, the same dispersion relation f l2"U]) for Alfven waves is obtained. However, the only 
other solutions present in the parallel direction are 



co = ±k^(3 + r). (26) 

These waves have frequencies which are purely real and correspond to undamped sound 
waves of double adiabatic model with isothermal electrons. 

For completeness, considering perpendicular propagation, it can be shown that the heat 
fluxes vanish and the solutions of the Landau fluid model with isothermal electrons are 
undamped magnetosonic waves with the dispersion relation expressed as 



± J 1 + ,„( ap + I) + (^). (27) 
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